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Abstract. Characterization of medium-range order in amorphous materials and its 
relation to short-range order is discussed. A new topological approach is presented here 
to extract a hierarchical structure of amorphous materials, which is robust against small 
perturbations and allows us to distinguish it from periodic or random configurations. 
The method is called the persistence diagram (PD) and it introduces scales into many- 
body atomic structures in order to characterize the size and shape. We first illustrate 
how perfect crystalline and random structures are represented in the PDs. Then, the 
medium-range order in the amorphous silica is characterized by using the PD. The 
PD approach reduces the size of the data tremendously to much smaller geometrical 
summaries and has a huge potential to be applied to broader areas including complex 
molecular liquid, granular materials, and metallic glasses. 
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Glasses have become increasingly useful and popular in materials engineering and 
industry. Nevertheless, its microscopic structure is not as clearly understood as 
crystalline solids. This is mainly due to the lack of long-range order (LRO). Hence 
the local structures, instead of LRO, have been extensively studied to characterize 
amorphous structures. In particular, short-range order (SRO), which describes atomic 
configurations of nearest neighbors, was well understood both experimentally and 
theoretically. However, it becomes clear that the larger-scale structures beyond neighbor 
atoms so-called medium-range order (MRO) are much more important than SRO in 
glasses. The amorphous structures are essentially characterized by possible connection 
of SRO that can build up continuous atomic configurations throughout the materials 
without any periodicity [1, 2, 3, 4, 5, 6, 7, 8]. 

The intrinsic structural features of MRO can be detected in several ways. For 
example, the split of the second peak in the radial distribution function is a sign for 
amorphous metals, whereas the first sharp diffraction peak of the structure factor 
is a sign for covalent amorphous solids [8]. Although these signs are evidence for 
the existence of MRO, their geometric origins are not clearly understood like the 
“periodicity” of LRO or the “chemical or packing order” of SRO. This is because MRO 
is definitely generated by many-body configurations, and the 2-body distributions such 
as the radial distribution or the structure factor does not properly describe its geometry. 

To describe the many-body atomic structure appearing in amorphous solids, angle 
distributions and statistics of topological quantities have been widely used so far. The 
bond-angle and torsion-angle distributions extract the partial information of 3-body 
and 4-body configurations, respectively [1], Even though these variables describe the 
configurations beyond the nearest neighbor, the scope of the scale is restricted within 
m(< 4)-body configurations. 

On the other hand, the topological quantification of amorphous structures such as 
the ring statistics and the Voronoi polyhedral analyses have been used to characterize 
atomic configurations of covalent and metallic amorphous solids, respectively. [1, 2, 
9, 10, 11, 12, 13, 14, 15, 16, 17]. These methods do not restrict the number of 
atoms to be considered, and are useful to classify the variety of many-body atomic 
structure in certain situations. However, they require building some geometric models 
from atomic configurations based on artificial criteria such as a threshold of the bond 
length. Furthermore, they do not explicitly extract the length scale. Hence they are not 
suitable for studying multi-scale phenomena, which are supposed to be important to 
determine MRO. Systematic methodologies to study many-body atomic structure with 
metric information are highly desired. 

In recent years, persistent homology and its graphical representation, called the 
persistence diagram (PD), have been invented in the applied mathematics community 
as a systematic method to extract geometric properties embedded in point cloud data 
[18, 19]. By regarding the atomic configuration as a point cloud data, we can apply 
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the method for geometric analysis of materials. Significantly, this method can deal with 
many-body configurations, and also can provide two length-scales, birth and death scales, 
which characterize the size of the n-dimensional holes. Namely, topological properties 
with metric information can be derived. Furthermore, PDs can be computed quite 
efficiently [20, 21, 22] so that it can handle complex configurations with a huge amount of 
atoms obtained by molecular dynamics (MD) simulations. Based on these observations, 
the PD approach is expected to be an adequate tool to characterize MRO in a systematic 
way. 

In this paper, we propose a methodology based on PDs to characterize MRO in 
amorphous solids. We first provide a brief introduction of PDs mainly for readers not 
familiar with them. Then we show two examples to illustrate the geometric meaning of 
PDs by using FCC crystal and a uniform random configuration. Then, we calculated 
PDs for the atomic configurations of amorphous silica obtained by MD simulations. 
Our method can be applied to a wide class of amorphous solids and complex liquids 
with atomic configuration data. Here, we chose amorphous silica, since it is one of the 
standard amorphous solids, and hence, is an appropriate model material for comparing 
our method to the other conventional tools such as ring statistics. Moreover, the 
extensivity of the Betti numbers, the structural hierarchy of the geometric objects 
appearing in PDs, and the decomposition method into single atomic components are 
also explained. From these arguments, we elucidate the geometry of MRO in amorphous 
silica. As a consequence, we conclude that the PD analysis is an appropriate method to 
describe MRO, complementary to the existing tools. 

2. Persistence Diagram 

The input for computing persistence diagrams is a pair A = {Q, R) of an atomic 
configuration Q = (aq, x 2 , xn) and radius parameters R = (rq, rq,..., rqv) for an 
./V-atom system. Here, oq G M 3 and rq are the three-dimensional position and the 
input radius of the i-th atom, respectively. Instead of setting a threshold by an 
artificial criterion and giving bonds between atoms, we introduce a ball Bi{a) = 
{x G M 3 | ||x — xf \ < rq(a)} at 5q with radius rq(a) = \Ja + rf for each i-th atom 
parameterized by a. Then, we study persistent topological features in the union of balls 
B{a) = IJili Bi(a) by changing the variable a. Typically, for a increasing family of a, 
this is a set of inflating atomic balls with centered at aq. 

For each a, n-dimensional holes in B(a) such as clusters, rings, and cavities for 
n = 0,1, and 2, respectively, can be identified by homology [23]. Then, for each n- 
dimensional hole c&, we can uniquely assign values a = bk and a = dk at which ck first 
appears and disappears, respectively. These values bk and dk are called the birth and 
death scales of the hole cy. The collection 

D n (A) = {(&fc, dk) G M 2 | k = 1, 2,...} (1) 

of all these birth and death scales {bk, dk) is called the n-dimensional persistence diagram 
of A. Here, we remark that the minimum value a m j n of a can be negative, i.e., 
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ft mm = —min {r 2 , rf, r^}, and the dimension of a (therefore bk and c4) is length 

squared. 

The PD is a two-dimensional scatter plot, ft follows from bk < dk that the 
points in the PD appear above the diagonal. From the construction, the birth and 
death points ( bk,dk ) are determined by the m-body atomic configuration x it , x i2 ..., x im 
comprising of the k-th hole Ck . Namely, the birth and death scales are represented by 
bk = b}-{x lx , Xj 2 ,..., x lrn ) and dk = dk(x n , x l21 ..., x lrn ). Different from the existing many- 
body variables such as bond angle 9k = 9k{x n , Xi 2 ,x l3 ) represented by m = 3, or torsion 
angle f>k = x i2 , x is , x i4 ) represented by m = 4, the number m of atoms to be 

considered is not restricted for the birth and death scales. Therefore bk and dk are 
reduced variables of metric properties of c*,. We use these variables to describe MRO. 

Intuitively, the birth scale bk indicates the maximum neighboring distance in 
Xi x ,Xi 2 , ...,Xi m , because Ck appears when the largest bond (n = 1) or cap (n = 2) in 
Ck is created. On the other hand, the death scale dk indicates the size of Ck, because Ck 
disappears when it is covered up in the inflated atomic ball model (Jili Bfa). 

We also introduce the life scale dk — bk that represents the robustness of the hole Ck 
under the change of the variable ft. Then, the points far from the diagonal in the PD 
have long life scales, meaning that they persist in a wide range of ft. In contrast, the 
points close to the diagonal, i.e., small life scales, correspond to holes that are sensitive 
to ft and can be considered noise. 

In terms of the characterization of MRO, the life and death scales are important. 
The former distinguishes the proper geometric objects from the topological noise, 
whereas the latter represents the sizes of holes. 

Persistent diagrams can be computed efficiently using freely available open-source 
software [20, 21], For the geometric model of the union of balls structure, we use 
the alpha shapes module in CGAL [22], We remark that even though the molecular 
simulations for bulk systems use periodic boundary conditions, we work with the point 
cloud data as a set of points in M 3 as input to the PDs. Hence, some differences will 
appear when we count the holes across the boundary, although these effects are negligible 
in a sufficiently large system. 

3. Typical Examples 

Before analyzing the covalent glass structure, we introduce two typical examples. For 
simplicity, we only consider monatomic systems. For a monatomic system, the input to 
the PDs is given by the configuration of its atoms and a uniform input radius r for all 
the atoms. By the construction of D n fA), PDs calculated with a given radius r are the 
same as ones with r — 0 after the transformation ( bk , dk) —* (bk + r 2 , dk + r 2 ). Hence, 
we choose the input radius to be zero, and hence the radius with the scale a is given by 
r(ft) = yfa. 

The first example is a periodic configuration as a model of perfect crystal. In 
particular, we choose the FCC crystal. In the unit cell, there are 4 atomic sites at 
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(0, 0, 0), (l/\/2, l/\/2, 0), (l/\/2, 0, l/\/2) and (0,1/\/2,1 /y/2) and the three primitive 
lattice vectors are set to be (\/2, 0, 0), (0, \/2,0) and (0,0\/2). Therefore the number 
density is p = y/2. These lengths are chosen to make the distance between nearest 
neighbors unity. We prepare a cubic cluster of the FCC crystal in such a way that the 
number of unit cells along each edge is n^, the length of each side is L — \[2n^ 1 and 
hence, the number of atoms is N = 4 n\. 

The top panels in Fig. 1 show the PDs for the FCC. Because of the periodicity, 
only a few points, each with high multiplicity, are found in each PD. The point at 
(0,1/4) in D 0 represents the length to contact the nearest neighbor. This is because 
each inflated atomic ball contacts its nearest neighbor atoms at r/o;) = 1/2, i.e, a = 1/4. 
In the same way, the point at (1/4,1/3) in D i represents equilateral triangle rings, and 
the points (1/3, 3/8) and (1/3,1/2) in D 2 represent regular tetrahedral and octahedral 
cavities, respectively. Here, the open circles at the points close to (1/2,1/2) in D\ 
and D 2 represent rings that are regular squares or isosceles right triangles and cavities 
of isosceles right triangular pyramids, respectively. These holes are parts of a perfect 
octahedron. As the radius of atomic ball becomes large, these holes are instantly covered 
right after they appeared. 

The periodicity also leads to high multiplicity. By a simple geometric consideration, 
the multiplicities at (0,1/4) in D 0 , (1/4,1/3) in Z/, and (1/3, 3/8) and (1/3,1/2) in D 2 
are calculated as mN.Neighb = 4n/ — 1, tutu. = 20n| — 24n| + 6nt, + l, mTetra. = (2 ul — l) 3 
and mocta. = 4 (nr,— l) 3 respectively. The coefficients of nj represent how many clusters, 
rings, and cavities exist in the unit cell in the bulk system. Specifically, in the unit cell, 
there are 4 particles, 8 tetrahedra and 4 octahedra. Every ring in the unit cell is a 
regular triangle, and is shared by a tetrahedron and an octahedron. Therefore we only 
need to count the number of rings in the octahedra. Each octahedron has 8 triangles 
on the surface, and thus there are 32 triangle rings in the unit cell. However, in each 
octahedron, one ring can be expressed as a linear combination [18, 23] of the other 7. 
Similarly, in each tetrahedron, one ring can be expressed as a linear combination of the 
other 3. Thus, we subtract 4 + 8 from 32, to get a total of 20 linearly independent 
rings in each unit cell. Terms proportional to nf,nL and the constant come from the 
boundary effect. Therefore as the number of atoms becomes large, the multiplicities per 
N converge as is seen in Fig. 2. 

In the case of the perfect FCC crystal, the points on the diagonal are due to 
numerical artifacts, and therefore their multiplicities do not scale as n\. If a small 
perturbation such as thermal fluctuation is added to the crystal and resolves the 
degeneracy, the points with high multiplicities in the PDs will split and distribute with 
finite width. Then, the distorted octahedral cavities generate the secondary rings and 
cavities close to the diagonal and yield the high multiplicities scaled to nf , as will be 
shown later (See Fig. 5). 

The second example is a uniform random configuration scattered in a cubic box. 
This is a model for the ideal gas. In this case, Q is composed of N points sampled from 
the uniform random distribution on [0, L] 3 , and hence the number density of the system 
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Figure 1. Dq,Di and Di are shown in the left, center and right panels respectively. 
The top panels correspond to the cubic cluster of the perfect FCC crystal with 
riL = 10. The point at (0,1/4) in D 0 represents the contact to the nearest neighbor 
atoms, (1/4,1/3) in D\ represents equilateral triangle ring, (1/3,3/8) and (1/3,1/2) 
in D 2 represent tetrahedron and octahedron cavities respectively. The bottom panels 
correspond to a uniform random configuration composed of N = 27000 particles 
scattered in the cubic box [0, L} 3 . Each point in the PDs for the perfect FCC crystal 
has high multiplicity, whereas PDs for the uniform random configuration have broad 
distributions. 


is p = N/L 3 . The bottom panels in Fig. 1 correspond to the PDs for the system with 
N = 27000. The PDs show a broad 2-dimensional distribution except for D 0 whose 
birth scales are zero. 

3.1. fd n {a): Betti Number 

As seen in the two examples, both the broadness of distribution and the multiplicity in 
the PDs well describe atomic structures of the bulk system. For further study of the 
multiplicity, we define the Betti numbers /3 n (a) as 
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Figure 2. Multiplicities at (bk,dk) in D n per the number of particle N for perfect 
FCC crystal are plotted as a function of N. These values converge to their bulk values 
like m(n)/n 3 , where m{n) is a third-order polynomial. 


for a given a. In Fig. 3, (3 n (a) for N = 4000, 32000, and 1715000 with p = \J2 for the 
FCC crystal (top panel) and N = 1000, 27000, and 125000 with p — 1 for the uniform 
random configuration (bottom panel) are plotted. 

For the crystal, /3 n (a) is a piecewise constant function, whereas /3 n (a) for the 
uniform random configuration is a single peak function. The plateau values of f3 n for 
the FCC crystal are represented by the multiplicities in D n . Namely, the plateau values 
of ApA, the first and second plateau of /? 2 are equal to mN.Neighb, m Tri., ^Tetra. +?«Octa.> 
and mocta. respectively. Therefore, they are proportional to N for the asymptotically 
large N. It is noteworthy that /^(a) cannot separate the number of tetrahedra and 
octahedra even though D 2 can. In this sense, PDs have richer information than the 
Betti numbers. 

For a uniform random configuration, the curves of j3 n (a) seem to be scaled to N as 
seen in the insets of the bottom panel in Fig. 3. To illustrate the scaling, we compare 
the peak values /3 n * = rna x a /3 n (a) against N and confirm that the peak values /3 n * are 
asymptotically proportional to N (Fig. 4). Therefore, for given a, /3 n (a) is an extensive 
variable for the system. However, only in the case n — 0 and for large a, this extensively 
fails because ft 0 (a) = 1 and is independent of the system size. 
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Figure 3. Betti numbers f3 n {n = 0,1,2) are plotted against the variable a for the 
perfect FCC crystal with number density p = \/2 (top) and the uniform random 
distribution with fixed number density p = 1 (bottom) for several N. The solid, 
dashed, and dots lines correspond to the /?o, /?i and fa respectively. Insets: Betti 
numbers per N are plotted as functions of a. 


3.2. £ n (b,d): Normalized Distribution for D r 


Because the Betti numbers are asymptotically proportional to N for the examples above, 

it is natural to introduce the two-variable normalized distributions for D n (A) 
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Figure 4. The peak values of /3 n (a) for the uniform random configuration are 
plotted against the number of the atoms N. These values are proportional to N 
in asymptotically large N. 


for n = 1 and 2. Here, the factor p -4 / 3 is introduced to make dimensionless. Because 
is a distribution function whose argument is ( b,d ), its horizontal and vertical axes 
are b and d, different from the scatter plot D n whose axes are bk and dk- In the case of 
n = 0 for the monatomic system with input radius r, every birth scale takes the same 
value bk = —r. Hence, we define £o as a function of d by 

—2/3 

keDn 

For the FCC crystal, as was seen in Fig. 2, because the multiplicity for each point 
asymptotically scales to N, converges as N become large. We also validate this 
convergence for non-zero temperature crystal obtained by MD simulation for Lennard- 
Jones (LJ) particles’ system. The parameters of the potential u(r) = 4e((<r/r) 12 —(a/r) 6 ) 
are set to be e = a = 1. Starting from the perfect crystal configuration, NPT (Nose- 
Hoover-Anderson) simulation at temperature T — 0.1 and pressure P = 1 in LJ unit 
has been carried out to obtain the equilibrium configuration of the FCC crystal. As 
a reference, we have also calculated PDs for the perfect crystalline configuration at 
T = 0, P = 1 obtained by energy minimization. According to the top panel of Fig. 5, 
is distributed with non-zero width around the spike corresponding to the perfect crystal. 
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Similar to the perfect crystal, the distributions for the £1 and £2 surfaces converge 
asymptotically large N. In the PDs, each geometric object forms one isolated island 
domain over the diagonal, which is characteristic for the crystalline solid. In addition, 
peaks close to the diagonal in £1 and £2 also converge asymptotically large N, different 
from the perfect crystal, therefore they are no more numerical artifact. Both rings and 
cavities in these regions are generated by the octahedra. This is a typical example, where 
one primary hole generates other secondary holes. In this case, the octahedral cavity 
is the primary hole and the rings and cavities close to the diagonal are the secondary 
holes. 

For the uniform random configuration, £, n (p 2 ^ 3 b, p 2 ^ 3 d) is a scaling functions. 
Intuitively, the birth and death scales become smaller as the number density p becomes 
larger. This effect is cancelled by replacing the arguments of £ n by (p 2 ^ 3 b, p 2 ^ 3 d), as 
seen in left panels of Fig. 6 for N = 6400, L — 40 and 10. The left panels in Fig. 6 
also show the statistical convergence for p = 1 with N = 1000, 64000, and 125000. For 
small N, only the points close to the diagonal appear because holes with short life scales 
are the majority in the random configuration. In contrast, a large N is necessary to 
detect holes with long life scales. The right panels in Fig. 6 show the scaling function for 
N = 125000. The derivation of this scaling function is explained in Appendix Appendix 
A. 


To grasp the meaning of the topological description using (2) and (3), we compare 
them to the coordination number ticn and the radial distribution function g(r). Recall 
that ncN and g(r) satisfy the following relation: 


ncN = 


drAnr 2 pg(r). 


(5) 


In Fig. 7, we see that £0 and 1 — Po/N show similar behaviors to g(r) and ticn- In fact, 
Po and £ 0 satisfy 

Po(oi) 


N 


= P 2/3 / dd£o(d), 


and by using Po(a m \ n ) = N, we obtain 


Po(a) 

N 



( 6 ) 


Different from the radial distribution g(r), £o(d) is determined by many-body atomic 
configurations. Actually, as long as the dimers are counted in the small scale, they 
yields the same distribution. However, as the scale becomes larger, such as the cluster 
composed of more than 2 atoms, these distributions deviate each other. For n = 1,2, 
the relation between £ n and P n is given by 

HP = y / 3 P d d [ a dbue d). (7) 

J OL J— OO 


Even though we show the extensivity of Betti numbers p n and the normalization 
of £ n only for the monatomic system, these properties are expected to be satisfied even 
for the multi-component system as long as the system is macroscopically uniform. 
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Figure 5. for Lennard-Jones FCC crystal of temperature T = 0.1 and pressure 
P = 1 in LJ unit. The convergence of for til = 50,20,5 (top panel). The contour 
map of £„ for n l = 50 (bottom panel). The spikes in top panel and the circle 
and triangles in bottom panel correspond to the configuration obtained by energy 
minimization, and therefore correspond to the configuration of perfect FCC crystal. 
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Figure 6. The convergence and collapse of for the uniform random configuration 
(left panels). For fixed number density p = N/L 3 = 1, the convergence of £„ is observed 
as N becomes larger (red, green and magenta surfaces). For fixed N = 64000, 
for L = 10 and 40 (green and blue surfaces) collapse as a function of (p 2 / 3 b, p 2 i 3 d). 
The contour maps of £„ (right panels) for the uniform random configuration with 
N = 125000, L = 50 as a function of (p 2 / 3 6, p 2 ^ 3 d). 

4 . Multi-Component System: Silica Glass 

In this section, we discuss how to apply the PD analysis to multi-component systems, 
and explain how to get multi-scale geometric information from the PDs. We choose silica 
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Figure 7. Correspondence of the relation between radial distribution function (right 
bottom) and coordination number (right top) and between £o (left bottom) and 
1 ~ fto/N (left top) for the FCC crystal of temperature T = 0.1 and pressure P = 1 in 
LJ unit. The length scale in £ 0 and 1 — /3 0 /N is rescaled to 2a 1 / 2 , which corresponds 
to the distance between two atoms in a 2-body distribution. 


glasses as examples. The atomic configurations of silica glasses have been obtained by 
the cooling MD simulation using BKS potentials [27, 25, 26]. 

4-1. Single-Component Analysis 

In a single-component analysis, the input to the PDs is given by the individual 
configurations of oxygen and silicon atoms extracted from SiCD atoms. Fig. 8 shows 
the PDs of the oxygen and silicon configurations, respectively. Similar to the previous 
section, we use the uniform input radius r = 0 for all atoms by using the invariance 
under the transformation. As we see later, this invariance property can be used to 
determine the atomic compositions of the holes in the multi-component analysis. 

We observed that both the birth and death scales of the silicon atoms are larger 
than those of the oxygen atoms. This is because the silicon atoms are distributed more 
sparsely than the oxygen atoms. Furthermore, we found several characteristic domains 
in the PDs: straight lines parallel to the death axis (fixed birth scale) in D\ , a curve 
departing from the diagonal in Tfi, and a spot domain I T around (b,d) = (2.4,2.6) in 
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Figure 8. PDs for the silicon (left) and oxygen (right) configurations in the amorphous 
silica. 

D 2 for the oxygen configuration. 

Comparing with the PDs for the uniform random distribution, these characteristic 
domains also encode geometric structures in the single-component configurations. For 
example, the vertical line determined by the fixed birth scale (b = bf) in D 1 show 
the existence of a geometric structure in the rings on the line. Recall that from the 
definition, the birth scale measures the distance to the neighboring atoms in the ring. 
This means that the rings on the straight line possess the fixed distance 2 r(a) = 2 y/R, 
to the neighboring atoms. On the other hand, the diversity of the straight line parallel 
to the death direction means a variety in the sizes of the rings. Moreover, we found that 
a ring on the curve close to the diagonal in D\ is generated by the ring on the vertical 
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Figure 9. £1 (left panel) and £2 (right panel) for an amorphous silica. Dq, D 1 and D 2 
are described in the inset. The characteristic curves are observed in £ 1 . 


line, which is similar to octahedron in FCC. We can also show that the spot domain It 
represents tetrahedra composed of 4 oxygen atoms each. 

4-2. Multi-Component Analysis 

4-2.1. Existence of Characteristic Curves Different from the single-component 
analysis, the PDs in this case depend on the choice of input radii parameter R = 
(ri,r 2 ,... ,tat). One of the ways to set R is using partial radial distribution functions 
of the atomic configuration. Namely, we investigate the first peak positions of pairwise 
distributions, and set the initial radius for each type of atom by solving linear constraints 
determined by these peak positions. For example, in the atomic configuration of the 
amorphous silica, we chose the smallest two peaks dsiO = 1-65 A and d 00 = 2.55A (Si-0 
and 0-0 pairs, respectively). Then, solving the two constraints rsi + ro = dsio and 
2ro = doo, we obtain the input radii rib = 0.375A and Vq = 1.275A. Fig. 9 shows the 
PDs of the amorphous silica for these input radii parameters. 

In the previous work [27], the PDs of silica glasses and the hierarchical geometric 
ring structures are studied in detail. We summarize here some of the results which are 
relevant to the later discussion in this paper. First of all, three characteristic curves, 
Op, Op and Co, and one band region B 0 are found in D l of Fig. 9. The rings on Op 
have the property The rings on Op have the property that they generate secondary 
rings when the rings get pinched as the parameter a increases [27]. Here P is named 
after primary and the mechanism is similar to the octahedron cavity in the FCC. The 
resulting secondary rings are located on Op, C-o, and Bq- 
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We also note that the cavity distribution in Dj for the amorphous silica shows no 
characteristic regions. In particular, the spot domain Jx disappears. This is because 
the silicon atoms are placed in the interior of the 4-oxygen tetrahedra. 

From the PDs, we can conclude that the ring structures are more essential than 
the cavities in the amorphous silica. This is consistent with the fact that existing 
methods such as ring statistics well characterize the amorphous silica. Hereinafter, we 
only discuss D\. 

4-2.2. Classification of Curves By changing the radius parameter R , we can obtain 
further information such as the composition of the rings. In Fig. 10, several D\ are 
shown for the different input radii. The radii are set to be rsi = rgj —0.2A, rg ; ,rg ; -|-0.2A, 
and rgj + 0.4A for the silicon, and tq = rf, — 0.4A, rfi — 0.2A, and rf, for the oxygen. 

As the input radius of oxygen becomes larger, the death and birth scale become 
smaller for fixed input radius of silicon r Si (the top left panel in Fig. 10). However, 
two regions Co and Bo overlap for the different choice of rsi after the transformation 
(i bk,dk ) —» (2 y/bj. + r'o,2^/dk + Tq) (the bottom left panel in Fig. 10). As we discussed 
in the single-component analysis (see the right panel in Fig. 8.), (bk + rQ,dk + r‘o) 
represents the transformation induced by the input radius Vq. Hence, the overlap on 
the transformed coordinates means that the birth and death scales for the rings on Co 
and Bo are determined by oxygen. 

Similarly, we found that C-p is invariant to the input radii after the transformation 
( b k ,d k ) ->■ (\Jbk + rgj + a Jb k + rQ,2y/d k + Tq) (the center panels in Fig. 10.). This 
means the birth scale is determined by the length between silicon and oxygen, and 
death scale is determined by the length between oxygen atoms. For Ct, we found the 
invariance via ( b k ,d k ) —> (\/b k + r'g ; + ^Jb k + Vq, y/b k + r§j + y/dfi + Tq) (the bottom 
right panel in Fig. 10), meaning that the both scales are controlled by the silicon and 
oxygen atoms. 

Using this procedure, we can classify the regions with respect to the atoms 
determining the birth and death scales. Furthermore, we can also investigate the 
compositions of the rings by hireling optimal rings [24], For example, it can be shown 
that the rings on Co and Bo consist of only oxygen atoms. 

4-2.3. Geometric Constraints Encoded in Curves In contrast to the 2-dimensionally 
broad distribution of the PDs for the random configuration, the curves indicate that 
there exist geometric relations in the atomic configuration. Concretely, the sharp 
distribution normal to the curve represents a geometric constraint on the rings, whereas 
the broad distribution along the tangential direction represents a variety of rings. 

As shown in Fig. 9, Cp lies parallel to the death axis on a foxed birth scale. It 
follows from the invariance with respect to the input radii after the transformation that 
the birth scale is determined by the bond between silicon and oxygen. Thus, we can 
conclude that the geometric constraint in Cp means there is a sharp distribution of the 
Si-0 bond length. On the other hand, the diversity of the arrangements of the rings 
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Figure 10. The PDs for several input radii are described. For fixed rsi = rg; = 0.375, 
ro = Tq — 0.2A,ro, and Tq + 0.2 A are described in the top left panel. For fixed 
ro = Tq = 1.275, rsi = *"so r Si + 0.2A and rg ; + 0.4A are described in the top 
right panel. The overlap of Cp is observed under the transformation (&*,, dp —> 
{y/bk + J'gi + \Jbk + Tq, 2yJdk + t-q) in the middle panels. The collapse of Co and 
Bo region is observed under the transformation ( bk,dk ) —t (2 yjbk + Tq, 2yJ dk + Tq) 
in the bottom left panel. The collapse of Cp is observed under the transformation 
0 bk,dk ) -t (\/bk +rg; + sjbk + ?o, yj d k + + yj d k +^o) in the bottom right panel. 
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Figure 11. Comparison between the ring statistics (bottom) and the distribution of 
death scales in Cp (top). Both of them represent the size distribution of the rings 
constructed by Si-0 network. 


is mainly controlled by the arrangement of oxygen atoms, and this is represented by 
the diversity along the death axis. The death scales in Cp describe the size of ring 
constructed by Si-0 network, and hence its diversity corresponds to the ring statistics 
[1, 8]. Therefore, as is shown in Fig. 11, the ring statistics for the shortest path rings 
[13] and the death scales in Cp show similar distribution. Even though both of them 
represent the size distribution of the rings, only the death scale encodes the metric 
information. 

Because the birth and death scales of Cp are controlled by the bond lengths between 
silicon and oxygen, Cp represents short range order. Fig. 12 shows both Cp and the 
Di for the distorted tetrahedra. Here, the distorted tetrahedra are constructed in such 
a way that tetrahedra in the local arrangement of SiC >4 require Si-0 and 0-0 bond 
lengths to be within the Erst peak of the distributions. Furthermore, their O-Si-O angles 
and Si-000 spherical angles are restricted within the 2.5 and 1.5 times the width of 
the deviation around their mean value, respectively, determined from the configuration 
obtained by the MD simulation. We found that the curve for the tetrahedron deviates 
from the Cp if we remove one of the restrictions above. From this observation, the 
agreement of the two PDs indicates that the geometric structure represented by the 
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Figure 12. Ct (red) and D i for the distorted tetrahedra (green) are described. 


curve Ct is a distortion of the SiC >4 tetrahedra. 

According to Fig. 9, the birth scales of Co and Bq are larger than the death scales 
of Ct- This means that they appear after all triangles in SiCb tetrahedra have been 
covered by the inflated atomic balls. In addition, recall that the rings on Co and Bq 
consist of the oxygen atoms in the primary rings at Cp. Thus, they are smaller than Cp, 
representing the order related to the neighboring of tetrahedra, the so-called Short Range 
MRO (SRMRO). The geometric constraint corresponding to Co shows the relationship 
between 0-0 length distribution and 0-0-0 angle distribution. For further details, we 
refer the reader to the paper [27]. 

Because PDs can encode many-body atomic structure, we can get a birds-eye view 
of the geometric hierarchical structure among the various kinds of SRO and MRO. The 
key point is that a curve in a PD encodes an order in the atomic configuration. Namely, 
the curve Cp represents SRO, whereas the curves Cp, Co and Bo show the geometric 
structure corresponding to MRO in amorphous silica consistent with existing methods. 

From the invariance to the input radii, we can convert the death scale c4 of each 
ring to the real space length 2 \J + Uq, which shows the diameter of the ring. Since 
the death scale for the MRO rings (Cp, Co and Bo) is determined by oxygen, this 
diameter is independent of the choice of the input radii. The paper [27] shows that this 
length scale reproduces the typical length scale of MRO corresponding to the first sharp 
diffraction peak. 
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5. Concluding Remarks 

Persistent homology is a suitable tool to extract many-body atomic structures, and 
successfully describes geometric features of MRO. The topological approach dramatically 
reduces the degree of freedom to represent the many-body structures, while keeping 
metric information and qualitatively important geometric features in MRO. 

Several concepts in PD for atomic structures are introduced by choosing two 
extreme cases: crystalline and random structures. For the crystalline solid, the periodic 
structure yields a few island supports in f n with high multiplicity, and the diagonal 
region corresponds to the secondary holes that represent distortion from the primary 
holes. The shape fluctuation for each hole is represented as a broadness of the peaks in 
each island. For the random configuration, there are no isolated supports in and a 
single broad support is found. These observations illustrate the power of PDs to express 
many-body atomic structures. 

As an example of network forming glass materials, £1 for the silica glass shows 
neither island support nor single broad support but rather characteristic curves that 
imply the existence of MRO. The primary curve in 0 describes the size distribution 
of rings, which is similar to the ring statistics. Combining with other curves in £i, 
geometric constraint in the MRO is encoded in the PD. With the aid of length scales 
encoded in PDs, a hierarchical and multi-scale many-body atomic structures are more 
successfully represented in an integrated manner compared to the existing methods. 
We believe that this method also has a great potential to describe the other disordered 
systems, such as complex molecular liquids, packed granular materials, and metallic 
glasses. 

The analysis by the persistent homology is an integrated method incorporating 
the various existing topological tools. The variables /3 0 and £ 0 are the many-body 
counterparts of the radial distribution function and the coordination number. The 
variable in Cp corresponds to the ring statistics. The variable £2 might correspond to 
the statistics of the Voronoi index, although not mentioned in this paper. Additionally 
Betti numbers fd n shows extensivity at least for monatomic systems. It is expected 
that the extensivity is still satisfied for the multi-component system. By using this 
behavior, the phase transition between macroscopically isotropic phases such as liquid- 
liquid transition might hopefully be more clearly categorized in the phase diagram. 


Appendix A. The Scaling of f n for Uniform Random Configurations 


Suppose that N points are sampled from the uniform random distribution in the 
D-dimensional box [0 ,L] D . Let the A-point configuration be denoted as = 

(xi,X 2 , ...,xn). Here x % is a D-dimensional position, and N is sufficiently large. Then, 
the distribution of the n-th persistence diagram D n (Q^) is given as 



= j2^^b k )s(d-d k ). 

k 


1 —'n 


(A.l) 
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Here, the (bk, dk ) are determined by the PD of the union of balls whose radius is 
ri(a) = \foi. The control parameter in the system is the number density p = N/L D . 

Set m an integer with m > 2. Because of uniform configuration, v/ ”* L shows 

“similar” behavior as m samples of Q^d f° r sufficiently large N, where we can neglect 
any boundary effect. Therefore, the distribution S n is simply multiplied by m 




(A.2) 


We expect that Eq. (A.2) be satisfied even when m is real number. The scaling in 
Eq. (A.2) leads to 


H„ (b,d:Q^)=Nf(b,d) + o[N). 


(A.3) 


Here f(b, d ) is the mean value of H n /N with respect to the uniform random distribution 
and is independent of N and L , for the sufficiently large N due to the central limit 
theorem. 

In particular, for the uniform random configuration, the birth and death scales 
change as (bk, dk) —>■ (oi 2 3 bk, oi 2 dk) under the transformation, Q^ nd —>■ \ = 

(atxi, ax 2 ,..., oixn)- In addition, aQ^d I s identical to Q^ nd L ■ Therefore 
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(A.4) 


By using Eqs. (A.2) and (A.4), we obtain 

, , n N,L\ v 4 5 6 7 8 9 10 _ ( b 

, b > d\ Q rnd 
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for the arbitrary u,v. Then, by setting u = 1/N, v = 1/L, 

>.4;QZt) = Wp t/D S„ ( p 2/D b,p^ D d-,Q") + o(N) 


(A.5) 


(A.6) 


is obtained. Here, the argument is formal notation. It should be rewritten by using 
Eq. (A.2) as 


H„ (6, d; Q^ 1 ) = Np 1 ' 0 / {p 2 ' D b, p 2 ' D d) + o(N), (A.7) 


where / = p 4//D f and / is dimensionless function. The bottom panels in Fig. 6 shows 
the collapse of f(p 2 / D b, p 2 / D d ) for D — 3. 
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